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Abstract 

A numerical study of confined jets in a cylindrical 
duct is carried out to examine the performance of two 
recently proposed turbulence models: an RNG-based 
K-e model and a realizable Reynolds stress algebraic 
equation model. The former is of the same form as 
the standard K-e model but has different model coef- 
ficients. The latter uses an explicit quadratic stress- 
strain relationship to model the turbulent stresses and 
is capable of ensuring the positivity of each turbulent 
normal stress. The flow considered involves recircula- 
tion with unfixed separation and reattachment points 
and severe adverse pressure gradients, thereby provid- 
ing a valuable test of the predictive capability of the 
models for complex flows. Calculations are performed 
with a finite-volume procedure. Numerical credibility of 
the solutions is ensured by using second-order accurate 
differencing schemes and sufficiently fine grids. Calcu- 
lations with the standard K-e model are also made for 
comparison. Detailed comparisons with experiments 
show that the realizable Reynolds stress algebraic equa- 
tion model consistently works better than does the stan- 
dard K-e model in capturing the essential flow features, 
while the RNG-based K-e model does not seem to give 
improvements over the standard K-e model under the 
flow conditions considered. 

1. Introduction 

The flow configuration considered in this paper is 
sketched in Fig.l. It involves an inner high speed round 
jet and a slowly moving annular stream, both interact- 
ing with each other. Because of turbulent entrainment, 
the jet increases its mass flux while spreading. This 
must be balanced by an equal decrease in the mass flux 
of the ambient flow. The decrease in the ambient ve- 
locity thus sets up an adverse pressure gradient which 

affects, in turn, the evolution of the flow. Depending 

on the ratio of jet to ambient velocities at the entrance, 

two different flow regimes occur in the downstream re- 

gion: if the ratio is small, the jet cannot consume all the 

ambient flow before reaching the duct wall so that the 

flow remains unseparated; if the ratio is large, the op- 

posite happens and further entrainment must create re- 

verse flow to maintain the total mass flux conservation. 


Further downstream, the flow completely loses its jet 
characteristics and degenerates eventually to the fully 
developed regime if the duct is long enough. These flow 
features can be found in many engineering apparatuses 
involving two flows of differing velocities, in particular, 
in combustion chambers and ejectors. Therefore, the 
understanding of confined jet flows is of great interest 
in engineering applications. 

From a turbulence modeling point of view, the con- 
fined jet flow also constitutes a valuable test due to 
its complicated flow features. It is noted that the flow 
past a backward-facing step is a standard test problem 
to benchmark the performance of turbulence models in 
complex flows. The confined jet flow has features simi- 
lar to those found in the backward-facing step flow such 
as recirculation with an unfixed reattachment point and 
severe adverse pressure gradient, and adds additional 
complexities arising from the motion of the separation 
point. 

Numerical calculations of confined jets have been 
reported by Gosman et al. 1 * * , Habib and Whitelaw 2,3 , 
Jones and Marquis 4 , Khalil et al. 5 and Zhu 6 . In 
these calculations, turbulence effects were represented 
either by the K-e model or by second-moment closures. 
However, the previous calculations were all made on 
very coarse grids and with the hybrid central /up wind 
scheme 7 that is highly diffusive in the presence of 
both convective dominance and flow-to-grid skewness. 
Therefore, they might largely be contaminated by nu- 
merical diffusion, and the results are far from conclu- 
sive. 

The purpose of the present study is to assess the 
performance of two recently developed turbulence mod- 
els in the confined jet flow. The models considered here 
are the RNG-based K-e model used by Speziale and 
Thangam 8 and the realizable Reynolds stress algebraic 
equation (RRSAE) model 9 , both within the framework 
of the two-equation formulation. The RNG-based K-e 
model is of the same form as the standard K-e model 10 
but assumes different model coefficients which are eval- 
uated by the theory. In the original version of the 
RNG K-e model, all the coefficients had constant values 
which have been shown by Speziale and Thangam 8 to 
be inappropriate. In the latest version of the RNG K-e 
model 8 , the model coefficient related to the production 
of dissipation term in the € equation is a function of 
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7 } y where 77 is the time scale ratio of the turbulent to 
mean strain rate. In the RRSAE model, the Reynolds 
stresses are calculated by a quadratic stress-strain re- 
lation. All the model coefficients in this relation are 
determined from the realizability analysis so that the 
model ensures the positiveness of the turbulent normal 
stresses. 

The test problem to be considered is taken from 
the experiment of Barchilon and Curtet 11,12 which pro- 
vides detailed experimental data. The flow can be 
characterized by the Craya-Curtet number C t which is 
the inverse square root of the total momentum, non- 
dimensionalized with the volume flux and the duct 
area 1 3. The experiment showed that recirculation oc- 
curs when C t < 0.96. Calculations are carried out with 
a conservative finite-volume method and on a numeri- 
cally accurate basis. As a common practice, the calcu- 
lation with the standard K-€ model is also included for 
comparison. Detailed comparisons with experiments at 
five Ct numbers clearly reveal the predictive capabilities 
of the models in these flows of great practical impor- 
tance. 


1 ) Standard K-€ (SKE) model 10 

^ = VtiUij + U jti ) - \K6 ih vt — Cp-j- (4) 

p O € 


l {r[UK ^ + ^ ]}+ ± {rl VK-(v^/-£ i} 


= G 


(5) 




= C lG i-C,^ (6) 

where G is the production term of the turbulent kinetic 
energy 


G— -[thI7i ) i-|-T22I^2 > 2+T’33^S,3 + ^ 12 (^ 1 ) 2 + ^ 2 ,i)] (7) 
p 

the velocity gradients Uij are calculated by 


2. Mathematical Formulation 


2.1 Governing Equations 

Incompressible, steady-state, turbulent flows are 
governed by the Reynolds-averaged continuity and 
Navier-Stokes equations. In the polar-cylindrical co- 
ordinate system (x, r) shown in Fig.l, the conservative 
form of these equations can be written as 


drU drV „ 
+ — = 0 


dx 


dr 


(i) 


d r trrrT dU d , /T ___ dU X1 

r dp , d dU t u 
p dx + dx ^ dx + p N 
5 . . dV ■ 

+ Tr 1 ' ^ + T )] 


( 2 ) 
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dr dr p r p 
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where U and V are the axial and radial velocities, re- 
spectively, p is the pressure, v is the kinematic viscos- 
ity and p is the density. The Reynolds stresses in 
Eqs.( 2) and ( 3) are calculated by using the following 
three turbulence models: 


dU dV V 

= ^,2 = IT, °m = -. 

dx dr r 


dU dV 

U\ 2 — , U 2 1 = 

1,2 dr y 2,1 dx 

and the model coefficients are 


C M = 0.09, Ci = 1.44, C 2 = 1.92, 
<tk = 1, = 1.3 


( 8 ) 


(9) 


2 ) RNG K-e model 8 . It is of the same form as the 
standard K-e model but uses the following coefficients 


C M 

Ci 


0.085, Ci = 1.42 - 


i»(l - vf 4.38) 


1 + 0.015I7 3 ’ 
1.68, <r K = <r t = 0.7179 


( 10 ) 


where 

KS 1 

V = — 5 = (2S tJ S,^/\ Sy = -{Uij + Uu) (11) 


3) RRSAE model® 


- = + U jti ) + Tij - \ KSij , i* = C„ — (12) 

P o € 


where 
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(13) 
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n = u k>l Ui, k , n = u klt u ki , (ie) 

* = n = (2n&n&) l/a . (16) 

% = {Ui,i-Vj,i)f 2 + 4c mii « m 

w m is the rotation rate of the reference frame, and the 
model constants are 

C T i = -4, C r2 = 13, C r3 = —2, A 2 = 1000. (17) 

In the work of Shih et al. 9 and Zhu and Shih 14 , the 
following two sets of values for A\ , 7 were tested 


Ai = 5.5, 7 = 0 

(18) 

Ai = 1.25, 7 = 0.9 

(19) 


and both of them have been found to give almost identi- 
cal predictions for the two backward-facing step flows. 
With Eq.( 19), the rotational effect of the mean flow 
enters into C M . However, we have found in the present 
work that the values in Eq.( 18) work better for the 
axisymmetric confined jets. Therefore the values in 
Eq.( 18) are taken here. The K and e in the RRSAE 
model are calculated with the same equations as in the 
standard K-€ model. 

2.2 Boundary Conditions 

Four types of boundaries are present in the calcu- 
lation; they are the inlet, outlet, axis of symmetry and 
solid wall. Among them, the inlet boundary conditions 
demand special attention because they have a consid- 
erable influence on the calculations 15 ’ 18 . Table 1 gives 
the inlet jet and ambient velocities taken from the ex- 
periment of Barchilon and Curtet 11 


Table 1. Inflow conditions 


c t 

Uj (cm/s) 

U a (cm/s) 

0.976 

1293.6 

84.48 

0.714 

1298.9 

60.72 

0.506 

1253.8 

39.81 

0.305 

1282.1 

21.86 

0.152 

1296.2 

7.42 


The Craya- Curtet number Ct is calculated by 

r _ Urn ( 20 ) 

~ m - Ui){d 0 /D o y + (U* - Ul)/2)V* K } 

where £> 0 =16cm, d 0 =1.2cm and U m is the mean veloc- 
ity of the section 

Ura = ( Uj - U a )(d 0 /D 0 ) 2 + U a (21) 

In the potential core and the ambient region (Fig.l), the 
velocity are uniform and the turbulence level is very low 


so that the flow may be treated as potential. However, 
between the potential core and the ambient region there 
exists a thin shear layer from which the turbulent en- 
trainment develops. The specification of the boundary 
conditions in this layer is nontrivial. In this work, the 
parabolic entrance region (PER) scheme of Zhu et al. 15 
is used. The PER scheme which allows the fine res- 
olution of the initial shear layer was developed on the 
assumption that although the flow as a whole is elliptic, 
there exists a short region near the entrance where the 
flow is parabolic. A parabolic calculation is first carried 
out over a short distance between x = 0 and x = x c , 
by using the following mixing length model 


_ dU 
uv = , 


I/J - C 2 (r 2 - n) 2 


dU 

dr 


( 22 ) 


where r\ and r 2 are the coordinates of the inner and 
outer edge of the initial shear layer (Fig.l) and C is an 
empirical coefficient given by 


C 2 = 0.0042 + 0.004 U a /Uj , 0 < U a /Uj < 0.2 (23) 


The results of the parabolic calculation are then used as 
the inlet conditions at x = x e for the elliptic calculation. 
The inlet values of K and e are calculated by 

K = -ut>/0.3 , e = 0.09if 2 /^ (24) 


It was found 6 that the PER scheme gives satisfactory 
predictions in the parabolic entrance region and the 
elliptic calculations were insensitive to i c provided that 
1 < x t jd 0 < 3. 

The outlet boundary is placed at x = 10D o where 
fully-developed flow conditions are assumed. Along the 
ft-ris of symmetry the normal velocity component and 
the normal gradients of the other variables are set to 
zero. The standard wall-function approach 10 is used to 
handle the wall boundary conditions. 

2.3 Numerical Procedure 

The transport equations ( 1), ( 2), ( 3), ( 5) and 
(6) can be written in the following generad form 

+ = (25) 

where ^ stands for U f V y K and e, and and 5^ 
are the corresponding diffusive coefficient and source 
term, respectively. For the momentum equations ( 2) 
and ( 3), 5^ also includes the cross-derivative diffusion 
terms and the quadratic terms in Eq.( 14). 

The system of equations ( 25) is solved with the 
finite- volume approach. It uses non- staggered grids 
with all the dependent variables being stored at the 
geometric center of each control volume. The velocity- 
pressure coupling is handled with the momentum in- 
terpolation procedure of Rhie and Chow 17 and the 
SIMPLEC algorithm of Van Doormal and Raithby 18 . 
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To ensure both accuracy and stability of numerical 
solutions, the hybrid linear/parabolic approximation 
(HLPA) scheme 19 is used to approximate the convec- 
tion terms of Eq.( 25). It has been shown 20 that the 
HLPA scheme of second-order accuracy works nearly 
as well as the third-order accurate SMART 31 and 
SHARP 32 schemes in eliminating the numerical diffu- 
sion while retaining the boundedness of numerical so- 
lutions. Considering a typical control volume centered 
at node C shown in Fig. 2. The HLPA scheme evaluates 
the value of <f> at the cell-face xv as follows 

K = U+<j>w + U~4>c + (26) 

where 




+ 


- 4>w) 


<i>w - 4>ww 
4>c - <t>ww 


Uw a w [<t>w - 4>c) 


4>c - 4>b 
<kw - 4 >e 


(27) 


Tr^)> u- = l-u: (U w * 0) (28) 


if | <j>c — 2 <f>w -j- <j>ww | < 1 4>c “ 4>ww | 
otherwise 

(29) 

if | <f>w - %<}>c + 4>e\ < I <f>w - 4>e\ / 3 qx 
otherwise ' ' 

It can be seen that Eq.( 26) is in fact the result of the 
first-order upwinding U+^w+U" <f>c with an additional 
term A <f> w added. The additional term may be viewed 
as an antidiffusive correction to the upwind scheme. 
The conventional central differencing scheme is used to 
approximate all other terms. The resulting discretized 
counterpart of Eq.( 25) can be cast into the following 
linearized form: 



A c <f>c = Ai<t>i + S, i = W 9 E, S, N (31) 

In formulating this equation, the convection terms cal- 
culated by the upwind scheme are coupled with the 
normal diffusion terms to form the main coefficients 
A*, while those calculated by Eq.( 27) are included in 
the source term 5. This way, the positivity of all the 
main coefficients is ensured so that the resulting coeffi- 
cient matrix will be always diagonally dominant. The 
system of equations ( 31) is solved with the strongly 
implicit solution algorithm of Stone 23 . The calculation 
results are considered converged when the maximum 
normalised residue of all the dependent variables is less 
than 0.5%. The details of the present numerical proce- 
dure are given in Rodi et al. 24 and Zhu 25 . 


3. Application 

All calculations were performed on the Cray YMP 
computer. The grid-dependency of solutions was first 
examined by using two convection schemes, HLPA 
and HYBRID (central/upwind differencing), and three 
grids consisting of 50x40 (grid 1), 86x50 (grid 2) and 
120x80 (grid 3) points, respectively. The HYBRID 
scheme that is highly diffusive in the presence of both 
convective dominance and flow-to-grid skewness has 
been used here mainly to highlight the importance of 
using higher-order accurate schemes. Test results ob- 
tained with the RRSAE model at C t =0.506 are shown 
in Fig. 3(a) for the axial velocity U- profiles, normalized 
by the mean velocity of the section U m , and in Fig. 3(b) 
for the turbulent shear stress tZr-profiles, both at the 
same downstream location x/D 0 = 1.875. It can be seen 
that the results of HLPA on the coarse grid 1 are al- 
ready very close to those on the fine grid 3 for both the 
U- and the uv-profiles, while significant differences ex- 
ist between the corresponding results of HYBRID. The 
HLPA results on the intermediate grid 2 can be consid- 
ered as grid-independent because the refinement from 
the grid 2 to the grid 3 produced differences too small to 
be seen on the graph. The HYBRID solutions, however, 
responded to the grid refinement in such a slow man- 
ner that they still had not reached the grid-independent 
stage on the finest grid. The numbers of iterations and 
CPU-time in minutes required for the calculations with 
HLPA were 196 and 0.2 on grid 1, 640 and 1.4 on grid 2 
and 1874 and 9.3 on grid 3. The calculations with HY- 
BRID took about 0.6~0.8 of these numbers. The grid 
2 and HLPA were used for all subsequent calcualtions. 

Fig.4 shows the variation of the centerline velocity 
U 0 with x and Cf. It clearly reveals the existence of 
the potential core characterized by the constant U 0 in 
the near-entrance region. Beyond the potential core, U 0 
decayed quickly, especially at small values of Ct. Both 
the SKE model and the RNG model predicted the same 
potential core length which was shorter than that pre- 
dicted with the RRSAE model at all the C t numbers. 
Since this length cannot be precisely determined from 
the first and second experimental points at each value 
of Ct, it is difficult to judge which model gives the bet- 
ter initial decay (x/D c < 1). For the ensuing decay, 
the RRSAE model gave the best agreement with the 
experiment while the RNG model produced large un- 
derpredictions. Figs. 5(a)- (c) show the axial mean ve- 
locity profiles at three C t numbers. All the three mod- 
els are seen to predict very well the upstream evolution 
of the flow. As for the downstream development, the 
results obtained with the RRSAE model remained in 
good agreement with experiments, while those obtained 
with the other two models deteriorated with the RN G 
model producing the largest discrepencies. The varia- 
tion of the ambient velocity XJ\ with x and C t is shown 
in Fig. 6. In the recirculation region, the ambient veloc- 
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ity has no physical meaning and is defined as the min- 
imum velocity (Fig.l) fox analytical convenience. The 
location where U\ is equal to zero corresponds to the 
separation or reattachment point. At C% =0.976, the 
calculated results are shown only up to xj D 0 — 1.875 be- 
cause the calculated 17-profiles have no uniform portion 
after this point (Fig.Sa). The calculated curves follow 
quite well with the experimental data upstream of the 
separation at all the C t numbers. The deviation occurs 
in the recirculation region. It should be pointed out 
that in the recirculation region, the computed velocity 
minimum* are all very close to the duct wall where the 
use of the wall-function as the boundary condition may 
constitute a source of error. It is also to be noted that 
Fig. 6 highlights considerably the difference between the 
computed and measured 77-profiles in the near- wall re- 
gion. This difference shown in Figs.5(b) and 5(c) is not 
as significant as in Fig. 6. Therefore, the RRSAE model 
result should be considered as satisfactory. 

The jet spreading can be characterized by the excess 
flow rate Qj and the effective width L They are defined 
by 


Q, = 2rf\u-U l )rdr, P = (32) 

Figs. 7 and 8 show the variation of Qj/Q and l/R with 
* and C u where Q and R are the total flow rate and 
the radius of the duct. As a result of the turbulent en- 
trainment, the excess flow rate increases first, reaches 
a maximum at the recirculation center where 77i has a 
minimum and then decreases. This variation becomes 
more marked as Ct decreases. Recirculation occurs 
when Qj is larger than Q . It can be seen from both 
figures that the calculations agree well with the experi- 
ments at larger C*, but the agreement deteriorates as C t 
decreases. It should be pointed out that the excess flow 
rate, due to its definition, is a quantity that is highly 
sensitive to the errors in the velocity profiles so that a 
small change in U\ 1 especially in the recirculation zone, 
will result in a large difference in Qj. Furthermore, the 
experimental uncertainty in the recirculation region in 
which the flow is highly perturbed is likely to be great- 
est. With due regard to these factors, the agreement 
between the calculations and experiments can be re- 
garded as reasonably good. Regarding the comparison 
among the three models, the RRSAE model again per- 
forms the best for both the excess flow rate and the 
effective width. 

Figs.9(a) and 9(b) show the predicted streamlines 
at the two values of C*=0.714 and 0.152. These fig- 
ures convey an overall view of the flow pattern. The 
upstream ambient flow was sucked in by the jet due to 
the turbulent entrainment. At C t =0.714, a small re- 
circulating bubble adhering to the duct wall occurred 
at the downstream location. When Ct was reduced to 
0.152, the recirculating bubble became very large, filling 


the most of the duct cross-section. The separation and 
reattachment points of the predicted recirculating bub- 
bles are compared with the experimental data in Fig. 10. 
The experiment indicated that as Ct decreased, the sep- 
aration point moved upstream while the reattachment 
point remained practically unchanged. The comparison 
shows that the RRSAE model gives the best predictions 
for both the separation and reattachment points. 

Fig. 11 shows the variation of the recirculating flow 
rate with * at C t =0.305 and 0.152. This is the inte- 
gral of negative velocities in each cross-section. The 
experiment indicated that the recirculating flow rate at 
Ct=0.152 is about 3 times larger than that at Ct=0.305. 
The results of the RRSAE model are in good agreement 
with the experiment while those of the standard K-€ 
model and the RNG model have substantial discrepan- 
cies. As for the maximum recirculating flow rate which 
is a critical parameter to characterize the performance 
of combustion chambers, the RRSAE model gave the 
same result as the experimental data at Ct =0.305 and 
a 9% overprediction at Ct =0.152 while the other two 
models produced larger overpredictions. It should be 
pointed out that results from different measurements 11 
for this quantity showed considerable scatter at small 
Ct numbers. The results of all the three models are 
within the experimental scatter. 

The variation of the pressure coefficient C v along 
the duct wall is shown in Fig. 12, where, C v is defined 
by 


„ _ &P-PUZ/2 

P ~ pUj / 2 

and Ap is the pressure difference relative to the en- 
trance. In the cylindrical duct, the evolution of the 
pressure is governed by the jet entrainment as well as 
the contraction and expansion of the flow caused by 
the recirculating bubble. The decrease in the ambient 
velocity (Fig.6) induced by the entrainment gives rise 
to an adverse pressure gradient, while the contraction 
of streamlines produces the opposite effect. These two 
mechanisms interact more intensely with each other as 
C t decreases, causing the pressure to vary little in the 
region upstream of the center of the recirculating bub- 
ble. In the downstream part of the recirculating bubble, 
the deceleration of the flow sets up an adverse pressure 
gradient the slope of which becomes steeper as Ct de- 
creases. Therefore, the ability to capture the location 
of the recirculation center will have a direct impact on 
the prediction of the pressure. By comparing Fig. 12 
and Fig.6, it can be seen that the three models cap- 
ture the steep pressure gradients in the same way as 
they capture the ambient velocity minimums. However, 
for the total pressure rise, an important parameter to 
the designer of jet pump devices, all the three models 
are seen to give the same results which are in excellent 
agreement with the measurement . 
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4. Conclusions 

A numerical study has been performed to assess 
two recently proposed turbulence models for confined 
jet flows. In order for the calculation to reflect the real 
performance of the models, an effort has been made to 
reduce numerical errors arising from the inlet boundary 
condition and numerical discretisation. The detailed 
comparison with the experiment definitively establishes 
the superiority of the RRSAE model over the standard 
K-e model in so far as the confined jet problem is con- 
cerned. However, this is not true for the RNG model 
at all the values of C t considered. 
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Fig. 4 Centerline velocity decay 
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Fig. 5 Axial mean velocity 
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(Notation as in Fig.4) 










Fig.6 Ambient velocity (Notation as in Fig. 4) 


Fig.7 Excess flow rate (Notation as in Fig.4) 
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